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Abstract. Transient diffusion equations arise in many branches of engineering and applied sci- 
ences (e.g., heat transfer and mass transfer), and are parabolic partial differential equations. It 
is well-known that, under certain assumptions on the input data, these equations satisfy impor- 
tant mathematical properties like maximum principles and the non- negative constraint, which have 
implications in mathematical modeling. However, existing numerical formulations for these types 
of equations do not, in general, satisfy maximum principles and the non-negative constraint. In 
this paper, we present a methodology for enforcing maximum principles and the non-negative con- 
straint for transient anisotropic diffusion equation. The method of horizontal lines (also known as 
the Rothe method) is applied in which the time is discretized first. This results in solving steady 
anisotropic diffusion equation with decay equation at every discrete time level, which is solved us- 
ing the methodology that has been recently proposed by Nagarajan and Nakshatrala (International 
Journal for Numerical Methods in Fluids, vol. 67, pp. 820-847, 2011). The proposed methodology 
for transient anisotropic diffusion equation will satisfy maximum principles and the non-negative 
constraint on general computational grids, and with no additional restrictions on the time step. We 
illustrate the performance and accuracy of the proposed formulation using representative numerical 
examples. We also perform numerical convergence of the proposed methodology. For comparison, 
we also present the results from the standard single-field semi-discrete formulation and the results 
from a popular software package, which all will violate maximum principles and the non-negative 
constraint. 



1. INTRODUCTION AND MOTIVATION 

Certain quantities (e.g., concentration of a chemical species and absolute temperature) naturally 
attain non- negative values. A violation of the non-negative constraint for these quantities will imply 
violation of some basic tenets of PhysicJ\ It is, therefore, imperative that such physical constraints 
are met by mathematical models and by their associated numerical formulations. In this paper, 
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1 There are some systems for which negative temperatures are possible (see Kittel and Kroemer [34]). Such cases 
are beyond the scope of this paper. 
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we shall focus on two popular transient mathematical models, in which physical restrictions like 
the non-negative constraint play a central role. The first model is based on Fick's assumption 
(commonly referred to as Fick's law) and balance of mass. Fick's assumption is a simple constitutive 
model to describe the diffusion of a chemical species in which the flux is proportional to the negative 
gradient of the concentration. The second model is based on Fourier's assumption and balance of 
energy, which describes heat conduction in a rigid conductor. Both these constitutive models 
combined with their corresponding balance laws give rise to transient diffusion-type equations, 
which are parabolic partial differential equations. 

There has been tremendous progress in Applied Mathematics for these type of equations with 
respect to existence and uniqueness results, qualitative behavior of solutions, estimates, and other 
mathematical properties \A9\ 120] . In particular, it has been shown that transient diffusion-type 
equations satisfy the so-called maximum principles [39]. It will be shown in a subsequent section 
that the non-negative constraint can be shown as a consequence of maximum principles under 
certain assumptions. Analytical solutions to several problems have been documented in various 
monographs (e.g., see references [1U|, 147] ) . However, it should be noted that most of these solu- 
tions are for isotropic and homogeneous media, and for simple geometries. For problems involving 
anisotropic and heterogeneous media, and complex geometries; finding analytical solutions is not 
possible, and one has to resort to numerical solutions. Obtaining physically meaningful numerical 
solutions for transient diffusion equation that satisfy maximum principles and the non-negative 
constraint is the main aim of this paper. It is well-known (and will be discussed in subsequent sec- 
tions) that many popular numerical schemes (including the ones that are based on the finite element 
method) do not satisfy maximum principles and the non-negative constraint. Even for isotropic 
diffusion, stringent restrictions on the time step and the computational mesh are necessary to meet 
these important mathematical properties. 

The usual approach of solving linear second-order parabolic partial differential equations under 
the finite element method is to employ Galerkin formalism for spatial discretization. Several theo- 
retical results (which include convergence proofs, a-priori estimates) for this approach can be found 
in the literature (e.g., see Reference [18]). But it has been adequately documented in the litera- 
ture that this approach will not satisfy maximum principles and the non-negative constraint (for 
example, see Reference [27], and also the discussion in Appendix). Thus, there is a need to develop 
new methodologies that will satisfy important mathematical properties like maximum principles 
and the non-negative constraint, and thereby improve the overall predictive capabilities of current 
numerical schemes. 

1.1. Maximum principles for diffusion- type equations in numerical setting. The first 
study on maximum principles in the context of finite elements can be traced back to the seminal 
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paper by Ciarlet and Raviart [16], which considered steady-state isotropic diffusion, low-order 
approximation, and simplicial elements. Since then, several other studies have addressed maximum 
principles for steady-state diffusion equation. A more detailed account of various works can be found 
in references [1SJ H3J [5U] . Although these papers have considered steady-state diffusion equation, the 
discussion in these papers is applicable to transient diffusion equations. A brief summary of these 
three papers is as follows. In Reference |45j . a non-negative methodology for mixed finite element 
formulation has been proposed for steady-state diffusion equation using techniques from convex 
quadratic programming. The paper also studied the effect of the non-negative methodology on the 
element local mass balance. In Reference [33], a methodology has been proposed for steady-state 
diffusion equation with decay that satisfies maximum principles and the non-negative constraint on 
general computational grids. (Note that the maximum principle for diffusion with decay is slightly 
different from the maximum principle with out decay.) This methodology will be utilized later in 
the present paper. In Reference [SD], a systematic study on the effect of high-order approximation 
on the violation of maximum principles and the non-negative constraint. In particular, it has 
been shown using numerical simulations that the violation of the non-negative constraint does not 
decrease with p-refinement. 

1.1.1. Maximum principles for transient systems. Several papers have also addressed maximum 
principles for transient systems (i.e., parabolic problems) in numerical setting. Herrera and Valocchi 
|28| have employed flow-oriented derivatives with backward Euler to obtain non-negative solutions in 
the context of finite difference and finite volume methods. One method that is commonly employed 
in the area of subsurface hydrology is by Chen and Thomee [13]. This method is based on the 
standard single-field formulation but employs lumped capacity matrix. (By the standard single- 
field formulation we refer to the formulation obtained by employing the semi-discrete approach using 
method of vertical lines at integral time steps, and Galerkin formalism for spatial discretization. 
See Appendix for more details of this formulation.) It is noteworthy that lumping capacity matrix 
approach is commonly considered as a variational crime [30j . Reference [7J also alters the capacity 
matrix to preserve positivity for parabolic problems but restricts to isotropic diffusion. Other 
notable works are [551 [53] [21] [19] , which all focused on getting restrictions on the mesh (and in 
some cases on the time step) to meet maximum principles. More importantly, they did not consider 
anisotropy, and such restrictions are not possible for anisotropic and heterogeneous medium. 

There are several papers that considered consistent capacity matrices, but derived restrictions on 
the time step to satisfy maximum principles [401 l60l [32] [271 [29] . A striking difference between the 
time step restrictions with respect to numerical stability and maximum principles is that numerical 
stability places an upper bound on the selection of the time step whereas maximum principles place 
a lower bound on the selection of the time step. The time step is selected based on the following 
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inequality: 




(1) 



where A^tabihty - g ^ e cr j|-j ca ^ ^ me s ^ e p t obtain stable results, and is the critical time step 



to satisfy maximum principles. It should be however mentioned that these works on deriving time 
step restrictions have considered one-dimensional problems or isotropic media, and these conditions 
are not applicable otherwise. To the best of our knowledge, none of the prior works presented a 
methodology for transient anisotropic diffusion equations to satisfy maximum principles and the 
non-negative constraint on general computational grids with no further restrictions on the time 
step. 

1.2. Our approach and main contributions of this paper. In this paper, we shall employ 
the Rothe method (or the method of horizontal lines) [57] to solve transient anisotropic diffusion 
equation. There are several papers in the literature that have employed Rothe method to solve 
parabolic equations [STJ [HJ [35j H] • These papers, except for Reference [27], did not apply the Rothe 
method in the context of maximum principles. Although Reference [27] addressed maximum princi- 
ples by using the Rothe method, but the formulation is restricted to isotropic diffusion. In addition, 
Reference [27] employed techniques from stabilized methods, which is different from the approach 
taken in this paper. In the proposed formulation, the temporal discretization using the Rothe 
method will give rise to inhomogeneous elliptic partial differential equation, which is solved using 
the approach presented in our earlier paper [33]. An attractive aspect of the proposed methodology 
is that there are no additional restrictions on the time step to meet maximum principles. 

1.3. An outline and notation used in this paper. The remainder of this paper is organized as 
follows. In Section [2j we present governing equations for transient anisotropic diffusion, and discuss 
maximum principles and the non- negative constraint. In Section [3l we derive a methodology for 
enforcing maximum principles and the non-negative constraint for transient anisotropic diffusion 
equation using the method of horizontal lines. In Section [U we illustrate the performance of the 
proposed formulation using representative numerical examples. Finally, conclusions are drawn in 
Section O with a discussion on plausible future works on enforcing maximum principles. 

The symbolic notation adopted in this paper is as follows. Repeated indices do not imply 
summation. (That is, we do not employ Einstien's summation convention.) We shall employ the 
standard notation for open, closed and half-open intervals [5]: 



(a, b) : = {x G R I a < x < b}, [a, b] := {x G R I a < x < b} 




(2) 
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Similar to our earlier paper [33], we shall make a distinction between vectors in the continuum and 
finite element settings. We also make a distinction between second-order tensors in the continuum 
setting versus matrices in the context of the finite element method. The continuum vectors are 
denoted by lower case boldface normal letters, and second-order tensors will be denoted by upper 
case boldface normal letters (for example, vector x and second-order tensor D). In the finite element 
context, we shall denote the vectors using lower case boldface italic letters, and the matrices are 
denoted using upper case boldface italic letters. For example, vector v and matrix K. Other 
notational conventions adopted in this paper are introduced as needed. 

2. GOVERNING EQUATIONS: TRANSIENT ANISOTROPIC DIFFUSION 

Let n C R nd be a bounded open set, where "nci" denotes the number of spatial dimensions. 
The boundary is denoted by d£l, which is assumed to be piecewise smooth. A spatial point is 
denoted by x G ft. The gradient and divergence with respect to x are denoted by grad[-] and div[-], 
respectively. Let t € [0,1] denote the time, where I > denotes the length of the time interval. 
The concentration of an inert chemical species is denoted by c(x, t). The (spatial) boundary is 
divided into two parts: T D and T N such that T D U T N = dft and r D n T N = 0. T D is that part of the 
boundary on which Dirichlet boundary condition (i.e., the concentration) is prescribed, and T N is 
the part of the boundary on which Neumann boundary condition (i.e., the flux) is prescribed. The 
unit outward normal to the boundary is denoted by n(x). The governing equations for transient 



anisotropic diffusion can be written as follows: 

^ ~ div[D(x)grad[c(x,i)]] = /(x,t) infix {0,1) (3a) 

c(x, t) = c p (x, t) on T D x (0, X) (3b) 

n(x) • D(x)grad[c(x,i)] = g„(x,t) on T N x (0,1) (3c) 

c (x,t = 0) = c (x) in ft (3d) 



where D(x) is the diffusivity tensor, /(x, t) is the volumetric source/sink, c p (x, t) is the prescribed 
concentration on the boundary, q p (x,t) is the prescribed flux on the boundary, and co(x) is the 
prescribed initial condition. The diffusivity tensor is symmetric, and is assumed to be bounded 
above and uniformly elliptic. That is, there exists two constants < £i < £2 < +00 such that 

£iy T y < y T D(x) y < ^y T y Vx e n and v y e R nd (4) 

The above initial boundary value problem given by equations (|3ap - (|3d|) is a linear parabolic partial 
differential equation. From the theory of partial differential equations, such equations are known 
to satisfy maximum principles under appropriate regularity assumptions on the input data and the 
domain [M1E2]- 
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Figure 1 . A pictorial description of parabolic cylinder J7j and parabolic boundary T%- 

Remark 2.1. It should be noted that a consequence of Fickian/ 'Fourier mathematical model is 
that a thermal/chemical disturbance at a point will be felt at other points instantaneously. This is 
because of the parabolic nature of the resulting partial differential equations. To put it differently, 
these mathematical models predict that the information travels at infinite speed, which is against 
the current accepted laws of Physics. Several modifications have been suggested in the area of heat 
conduction to have finite speeds for thermal disturbances, and most of these models are hyperbolic 
partial differential equations. Some notable works on this topic are by Maxwell [3H] ; Catteneo 
and Gurtin and Pipkin [24] . A more detailed discussion with respect to finite speed thermoelasticity 
can be found in Reference [31]. It is noteworthy that hyperbolic partial differential equations do 
not possess maximum principles "similar" to the ones possessed by elliptic and parabolic partial 
differential equations. This area of research is far from settled, and is beyond the scope of this 
paper. 

2.1. Maximum principles for parabolic equations. Maximum principles for parabolic partial 
differential equations can be traced back to Levi [36] and Picone [52] . A brief history and other 
references on maximum principles for parabolic partial differential equations can be found in the 
book by Protter and Weinberger [33]. Herein, we shall employ an approach similar to that of 
Nirenberg [46] . Before we state a maximum principle for linear parabolic partial differential equa- 
tions, we shall introduce relevant notation and definitions. The parabolic cylinder is defined as 
fiz := £1 x (0,1). The parabolic boundary is defined as follows: 

T% '■= |(x, t) E Cl% xGffiori=o} (5) 

The parabolic cylinder and parabolic boundary are pictorially described in Figure [TJ Let C m (Vt) 
denotes the set of functions defined on O that are continuously differentiable up to m-th order. We 
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shall introduce the following function space with differing smoothness in the x- and t- variables: 

C?(fiz) := |c : Ox ->• R | c, — , — € C(fi x ); t, j = 1, • • • , ndj (6) 

Theorem 2.2 (maximum principle) . Letc(x.,t) G C^(r2j)nC(Oi) satisfy dc/dt— div[D(x)grad[c]] > 
in Qj. TTien c(x, i) achieves its minimum on the parabolic boundary ofilx- That is, 

min c(x, i) = min c(x, i) (7) 
(x,t)Gf2i (x,t)er x 

Proof. A proof can be found in standard books on partial differential equations (e.g., see [Ml ESI 

[20]). □ 

Remark 2.3. The above maximum principle implies that if one has volumetric source everywhere 
and at all times (i.e., /(x, i) > 0) then the minimum will occur on the boundary of the domain 
or in the initial condition. A logically equivalent statement of the above theorem can be written 
as follows: If c(x, i) satisfies dc/dt — div[D(x)grad[c]] < ; the maximum occurs on the parabolic 
boundary. That is, 

max c(x, t) = max c(x, t) (8) 
(x,t)ef2z (x,t)er x 

Maximum principles play a central role in the study of partial differential equations. Many 
uniqueness theorems and powerful estimates for elliptic and parabolic partial differential equations 
utilize some form of maximum principles [23\ EE9] . Maximum principles also have important physical 
implications in mathematical modeling, as they place restrictions on physical quantities. One such 
implication is the non-negative constraint. We now show that, under certain assumptions, the 
non- negative constraint is a consequence of the maximum principle given by Theorem 12.21 For the 
present discussion, let us assume that T D = dtt (that is, we prescribe Dirichlet boundary conditions 
on the whole boundary). If f(x,t) > (i.e., we have volumetric source), c p (x, t) > (i.e., we have 
non- negative prescribed Dirichlet boundary conditions on the whole boundary), and co(x) > 
(i.e., we have non- negative prescribed initial concentration); then the maximum principle given by 
Theorem 12.21 implies that the quantity c(x, t) is non-negative in the whole domain and at all times. 
That is, 

c(x, t) > Vx G and Vt G [0, 1} (9) 

It should be noted that the above discussion on maximum principles and the non-negative constraint 
is in continuum setting. For most practical problems (which will involve complex geometries and 
spatially varying coefficients), it is not possible to find analytical solutions. Therefore, one has to 
resort to numerical solutions. This leads to the following questions, which are central to this paper. 
Whether numerical formulations satisfy maximum principles and the non-negative constraint for 
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transient diffusion equation. If so, under what conditions? If not, is it possible to fix a given 
numerical formulation to meet these important principles? This area of research is popularly 
referred to as discrete maximum principles. 

Remark 2.4. Some recent efforts \37\ |4"51 |4"3] have addressed similar questions with respect to 
maximum principles and the non-negative constraint, but all these studies have considered steady 
diffusion equation. 

2.2. Discrete maximum principles. The discrete analogy of maximum principles is commonly 
referred to as discrete maximum principles (DMP). Some main factors which affect numerical 
solutions with respect to discrete maximum principles are: 

(i) topology of the domain (e.g., shape of the domain, features like holes in domain), 

(ii) type of mesh (e.g., Delaunay, well-centered, structured vs. unstructured), 

(iii) element type (simplicial vs. non-simplicial elements), 

(iv) mesh size (i.e., aspect ratio), 

(v) medium properties (e.g., anisotropy, heterogeneity), 

(vi) order of approximation (i.e., low-order vs. high-order), and 

(vii) temporal discretization (e.g., time stepping scheme, selection of the time step). 

The first six factors are equally applicable to steady anisotropic diffusion equation. Systematic 
studies on the effect of first five factors on maximum principles and the non-negative constraint 
can be found in references [451 l4"3l I41j . Reference [50] discusses in detail about the sixth factor. 
The last factor (in combination with other six factors) is the subject matter of this paper. 

This leads to the problem statement of this paper: Develop a finite element methodology for linear 
transient tensorial diffusion equation that satisfies maximum principles and the non-negative con- 
straint on general computational grids for low-order finite elements with no additional restrictions 
on the time step. To the best of our knowledge, such a methodology does not exist in the literature. 
In the next section, we shall extend the optimization-based methodologies that are presented in 
references [45 \ |4"3] for steady diffusion equations to transient diffusion equation. We shall explicitly 
enforce constraints on the nodal concentrations to satisfy maximum principles and the non-negative. 
We shall restrict to low-order finite elements, which include two-node line element, three- node tri- 
angular element, four-node quadrilateral element, four- node tetrahedron element, eight-node brick 
element, and six-node wedge element. However, it should be noted that the proposed methodology 
is not applicable to high-order elements, as enforcing non-negative constraints at nodes does not 
imply non- negative concentrations throughout the domain for high-order elements (e.g., three-node 
line element, six-node triangular element) |5Uj . 
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3. PROPOSED METHODOLOGY: DERIVATION AND IMPLEMENTATION 

DETAILS 

Herein, we shall employ the method of horizontal lines (also known as the Rothe method) |57j as 
opposed to the commonly employed method of vertical lines [30] . The method of horizontal lines 
is a discretization sequence in which the time is discretized first followed by spatial discretization. 
To this end, we shall define two sets of time levels: integral and weighted time levels. The time 
interval of interest [0,2] is divided into N non-overlapping subintervals such that 



N 



[0,2] = \J[t n -x,t n ] 



(10) 



n=l 



where t n (n = 0, • • • , N) are referred to as integral time levels. For convenience, we shall assume 
that the time step At to be uniform, which implies that 



At = — and t n = nAt 



(11) 



However, it should be noted that the proposed methodology can be easily extended to non- uniform 
time steps. We shall apply the method of horizontal lines at weighted time levels, which are defined 
as follows: 



t n+v := (1 - rj)t n + nt n+1 



(12) 



where the parameter n £ [0, 1]. The concentration and its rate at integral time levels are respectively 
denoted as follows: 



C W(x) 



c(x,t = t n ) 
dc 



t;(»)(x) = -(x,t = t n ) 



(13a) 
(13b) 



The following notation is used to denote quantities at weighted time levels: 



C (n+V) 




:= (l-77) C W(x) 




i c(x, tn+r; ) 


(14a) 


V (n+V) 


;x) 


:= (1 -r])v (n \x) 


+ r?u (n+l)( x ) r 




(14b) 


c {n+n) 


x) 


• — Cp(x, t n _|_^) 






(14c) 




x) 


:= /( x > tn+r)) 






(14d) 






:= q p (x,t n+v ) 






(14e) 



3.1. Derivation. In designing the proposed methodology, attention will be exercised on two dif- 
ferent aspects. The first aspect is to make sure that the non-negative constraint and maximum 
principles are preserved after both temporal and spatial discretizations. The second aspect is to 
achieve numerical stability in solving the resulting differential-algebraic equations. As we shall see 
in subsection l3.21 we will be adding additional equations in the form of lower and upper bounds (i.e., 
inequality constraints). This implies that we will be dealing with differential-algebraic equations. 
It is important to note that numerical time integration schemes that are designed for ordinary 
differential equations may not be stable and accurate for solving differential- algebraic equations. 
This point has been discussed adequately in the literature (e.g., see references [H [251 [22] )• An 
important work on numerical time integration of differential-algebraic equations is by Petzold [51] , 
and the title of this paper ("Differential/algebraic equations are not ODEs") succinctly summarizes 
the above discussion. 

We shall employ the generalized-a method for temporal discretization. The generalized-a method 
was first proposed for second-order transient systems in Reference [15], and later modified for first- 
order transient systems in Reference [33] , After applying the generalized-a method to the governing 
equations (f3a|) - (f3c|) . we obtain the following equations: 

( x ) _ div[D(x)grad[c( n+Q / ) ]] = /("+<*/) (x) in ft (15a) 
c (»+°/)(x) =4 n+Q/) (x) onT D (15b) 

n(x)-D(x)grad[c (ri+Q / ) ] =4 n+Q/) (x) on T N (15c) 

where the parameters a m , ctj € [0,1]. In addition, we have the following relationship: 

c (n+l) ( x ) = c (n) ( x ) + At ((1 - 7 y«) (x) + 7 t;(™+ 1 ) (x)) (16) 
where the parameter 7 £ [0, 1]. The initial condition takes the following form: 

c (°)(x) = c (x) inO (17) 

Remark 3.1. Many popular time stepping schemes are special case of generalized-a method. For 
example, forward Euler (a m = l,a/ = 1,7 = 0), trapezoidal rule (a m = l,a/ = 1,7 = 1/2), and 
backward Euler (a m = l,af = 1,7 = 1). 

Herein, we shall take a m = 7. This selection is intended to inherit the non-negative property for 

the resulting time discrete equations. The time discrete equations in terms of concentration take 
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the following form: Find c ("+ Q /)(x) such that we have 

_L- c («+ Q /)(x) _ div[D(x)grad[c (n+a / ) l] = /( n +"/)(x) + — L-c< ft >(x) in fi (18a) 

a/At ay At 

c ( n+Q /)(x) = 4 n+Q/) (x) onr D (18b) 

n(x)-D(x)grad[c( n+Q /)] =^™ +a/) (x) on T N (18c) 

The above boundary value problem (|18ap - (|18cp is a second-order inhomogeneous elliptic partial 
differential equation with Dirichlet and Neumann boundary conditions. Specifically, equation (j!8a|) 
is the well-known steady-state anisotropic diffusion equation with decay, as a/ At will be always 
positive. The decay coefficient can be identified as l/(a/At), and the volumetric source term is 
/ ("+"/) (x) + A t c ( x ) • This boundary value problem is also known to satisfy maximum principles 
and the non-negative constraint. The selection a m = 7 made it possible to preserve maximum 
principles and the non-negative constraint by ensuring the decay coefficient to be positive, and the 
volumetric source at discrete time levels to be non-negative. It should be noted that any arbitrary 
temporal discretization will not preserve maximum principles and the non-negative constraint, 
which will be illustrated in Appendix. 

Remark 3.2. Diffusion equation with decay is of following form: 

a(x)c(x) - div[D(x)grad[c]] = /(x) (19) 

with a(x) > 0. If a(x) < 0, the equation is referred to as Helmholtz equation. It should be noted that 
Helmholtz equation does not have a maximum principle similar to the one possessed by diffusion 
equation with decay [23] . 

Recently, Nagarajan and Nakshatrala [43] have proposed a procedure for enforcing maximum 
principles and the non-negative constraint for steady diffusion with decay equation, which we shall 
modify to solve equations (|18a.[) (|18c|) . We start by applying Galerkin formalism to equations 
(|18aj) — (|18cp . The corresponding weak form takes the following form: Find c^ n+a f\x) £ V n + a f such 
that we have 

/ w(x)—^—c( n+a f\x) dn+ [ gradH • D(x)grad[c (n+a / ) ] dfi 
Jn OLjbt J n 

= [ w(x)/( n+Q /)(x) dn+ ! w{-k)-^— c^(x) dn 
Jn Jn a f At 

+ [ w(x)^ n+a/) (x) dn Vw(x)eQ (20) 
Jr N 
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where the function spaces V n + af and Q are defined as follows: 

V n+af := {c(x) G H 1 ^) | c(x) = 4 n+Q/) (x) on T D } (21a) 

Q := { w ( x ) e ff 1 ^) I tw(x) = 0onr D } (21b) 

After executing the usual steps of the finite element method, the above weak form (|20p can be 
converted to a system of linear equations of the following form: 

Kc (n+a/) = f(n+a f ) ^ 

where "ndofs" denotes the number of (free) degrees- of-freedom, c ( n+a /) £ R nd °f s denotes the 
unknown vector containing nodal concentrations at the weighted time level t n + af , f( n+a f) g ^ndofs 
is a known vector, and K is a symmetric and positive definite matrix. It will be shown in a 
subsequent section that the finite element solution obtained by solving the system of linear equations 
(|22|) may not satisfy maximum principles and the non-negative constraint. Using optimization- 
based techniques, we now modify the above solution procedure to meet these important physical 
constraints. 

3.2. Enforcing maximum principles and the non-negative constraint. We shall denote the 
standard inner product on finite dimensional Euclidean spaces by (•;•)• We shall use the symbols 
H and y to denote component- wise inequalities for vectors. That is, for given any two (finite 
dimensional) vectors a and b 

a <b means that < 6j Vi (23) 
Similarly, one can define the symbol y. The optimization problem can then be written as follows: 
minimize - ( C ( n+a /); Kc^ n+a A - (c^ n+a ^; f( n + a A (24a) 

subject to it^l^c^^cL^l (24b) 

where 1 is a vector containing ones of size ndofs x 1, and c m i i ^ Q/ ^ and Cm^' are respectively the 
lower and upper bounds. For enforcing maximum principles, c^^" 1 ^ and Cmi"^ can be taken as 
follows: 

c m^ af) '■= min \ mi S c o( x )> niin c{, n+a/) (x) \ (25a) 

(n+o/) J / \ (n+a f ), ,\ /nxu\ 

Cmax :=max<max cnlx , max c„ x > (25b 
For problems involving only the non- negative constraint, one can employ the following: 

= and ctt: s) = +oo (26) 
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Alternatively, for enforcing the non- negative constraint, one can replace the constraint (|24b|) with 
the following: 

-< c (n+a/) (27) 

where denotes the vector of size ndofs x 1 containing zeros. It should be noted that the above 
optimization problem (|24p belongs to quadratic programming. Since, for the problem at hand, the 
matrix K is positive definite (which makes the objective function (|24ap convex) the optimization 
problem belongs to convex quadratic programming [9]. 

Remark 3.3. It is important to note that solving a problem in quadratic programming, in gen- 
eral, is NP-hard [59j . This means that there is no efficient algorithm to solve a general quadratic 
programming optimization problem in polynomial time. However, a convex quadratic programming 
optimization problem can be solved in polynomial time, and several efficient algorithms are available 
in the literature [48 1 [6Tj [9] . Some popular packages that can handle convex quadratic programming 
optimization problems are MATLAB [2\, GAMS p], TAO 02], and DAKOTA [3]. 

One can then obtain the nodal concentrations at integral time levels as follows: 

c (n+l) = V JJ / 2g x 

a/ 

Although c( n+a /) y 0, the nodal concentrations at integral time levels based on equation (|28p need 
not be non- negative if a/ 7^ 0. To put it differently, one is assured of satisfying maximum principles 
and the non-negative constraint under the proposed methodology if a m = 76 (0, 1] and atf = 1. If 
needed, calculate nodal rate of concentrations using the following expression: 

r (n+l) _ Jn) _ (-1 _ ^)Af v (n) 

V ( n+ V = - - 1 7J (29) 

jAt K ' 

To obtain stable and accuracy results for the rates, one need to choose 7 > 1/2. The theoretical 

basis for this is given in Reference [H], in which it has been shown that for 7 < 1/2 the results at 

integral time steps t>( n+1 ) will not be bounded when calculated from the results at weighted values 

•y( n +7)_ In particular, see Proposition 5.2, and Figures 2 and 3 in Reference |44j . The various 

steps involved in the numerical implementation of the proposed methodology to satisfy maximum 

principles and the non-negative constraint are summarized in Algorithm [H which could serve as a 

quick reference during computer code design and implementation. 

4. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we shall illustrate the performance of the proposed methodology for enforcing 

maximum principles and the non-negative constraint using several canonical problems. We shall 

also perform numerical convergence studies on the proposed methodology. We shall restrict our 
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Algorithm 1 Implementation of the proposed methodology based on aif = 1. 



1: Input: Initial condition c(x), Dirichlet boundary conditions c p (x, t), Neumann boundary con- 
ditions q p (x,t), time step At, total time of interest I, a m = 7 £ (0, 1]. 

2: Construct initial nodal concentrations c^ ' 

3: Set cW < — c(°), t< — 0, n i — 

4: while t < X do 

5: Calculate and Cma^ (see equations 

6: Call non-negative solver to obtain c^ n+1 ^ 

minimize -( C ( n+1 >; Kc^+V) - ( c («+l) ; /(«+!)) 

c (n+l)g:J;™<2o/s 2 

subject to =< =< 



7: If needed, obtain rate of nodal concentrations at integral time levels (but need to choose 
7 > 1/2 to obtain stable results for the rates) 

{n+1) _ C ("+P - eg - (1 - TjAtejj 



7At 



8: Set c( n ) < — c( n+1 ),t^ — t + At, n< — n + 1 
9: end while 



numerical studies to one- and two-dimensional problems. It should be, however, noted that the 
proposed methodology is equally applicable for solving three-dimensional problems. We do not solve 
any three-dimensional problem here as, in comparison with one- and two-dimensional problems, 
there are no additional difficulties other than the usual book keeping that is associated with most 
three-dimensional problems. In all our numerical simulations we have employed low-order finite 
elements, and have taken a m = at = 1. It is assumed that 7 = 1, unless stated otherwise. 

4.1. One-dimensional problem with uniform initial condition. The following one-dimensional 
problem is taken from Reference |10] , which is also used as a test problem in Reference [27] in the 
context of discrete maximum principles. The computational domain is := (0, 1). The governing 
equations of the test problem take the following form: 

^-%^ = inO x := ( 0,l)x(0,X) (30a) 

=°i c(x = = Vt€(0,Z] (30b) 

CJX. 

c(x, 0) = 1 Vx G [0, 1] (30c) 
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The analytical solution can be written as follows: 



4 ~ (-1)™ 
c(x, t) = — > -. r exp 



(2n + l) 2 7r 2 t 



(2n + 1)ttx 
cos (31) 

n=0 

The analytical solution is bounded between zero and unity. In the numerical simulation, we have 
divided the computational domain into five equal linear finite elements, and have taken the time 
step to be At = 0.001 s (which is chosen arbitrarily). Figure [3] compares the analytical solution with 
the numerical solutions obtained using the single-field formulation and the proposed methodology 
for 7 = 1.0 and 7 = 0.1. The single-field formulation violates the maximum principle, as the 
obtained numerical solution is greater than unity. The proposed methodology satisfies the maximum 
principle, and gives stable results for concentrations under both 7 > 1/2 and < 7 < 1/2 cases. 
However, for stable and accurate results for rate of concentration, one needs to employ 7 > 1/2. 

Remark 4.1. For this problem, the initial condition is not compatible with the boundary conditions. 
That is, the (homogeneous) Dirichlet boundary condition at the right end of the domain at time 
t = is not equal to the initial condition. Hence, there is no classical solution to the initial 
boundary value problem given by equations (|30ap ~ (|30cp in the sense that c(x, t) G C 2 (Qx) nC(fix)- 
The analytical solution given in equation (|3ip should be interpreted in Lebesgue measurable sense. 

4.2. One-dimensional problem with non-uniform initial condition. Consider the following 
simple one-dimensional problem with homogeneous forcing function. This problem is a modification 
to one of the examples given in Reference [T7]. The initial boundary value problem can be written 
as follows: 

dc(x,t) d 2 c(x,t) „ , , . . 

— dt — tt^ 1 := ( ' ) x ( ' ^ ( } 

c (x = 0, t) = c(x = 1, t) = Vt G (0,2] (32b) 



c(x,0) 



1 if x G [a, b] 
otherwise 

The analytical solution to the above problem is given by 



(32c) 



2 00 1 

c(x, t) = — ~ (cos(mra) — cos(n-/r6)) sin(n-7rx) exp[— n 2 -K 2 t] (33) 

11 71=1 n 

In this paper, we have taken a = 0.4 and b = 0.6. 

Figure U shows that the numerical solution from the proposed methodology compares well point- 
wise with the analytical solution, and satisfies the maximum principle and the non-negative con- 
straint. In Figure [5J we have shown the numerical convergence of the proposed methodology with 

the standard single-field formulation in L2- norm an d -f/^-seminorm, which show convergence in 
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integral sense. The convergence study is carried out by employing simultaneous spatial and tem- 
poral refinements satisfying the condition At oc (Ax) 2 . The coarsest mesh has 100 nodes, and the 
corresponding time step used for this mesh is At = 0.05 s. 

FiguresE]and[7]show, respectively, the variation of the minimum concentration and the maximum 
concentration in the domain with respect to time under the standard single-field formulation. Note 
that for this problem the minimum concentration should be zero, and the maximum concentration 
should be unity. Clearly, the results from the standard single-field formulation violated both the 
upper and lower bounds. However, the extent of the violation decreased with time, which is 
expected as diffusion is a dissipative process. Figures [8] shows the effect of mesh refinement and 
the selection of small time steps on the discrete maximum principle for the standard single-field 
formulation. For a given mesh, the extent of the violation will be greater for smaller time steps. On 
the other hand, for a given time step, the extent of the violation decreases with mesh refinement, 
which is will not be the trend in the case of anisotropy. Figure E] shows the performance of the 
proposed methodology for two different time steps and for two different meshes. The proposed 
methodology satisfies the discrete maximum principle even on coarse meshes and for small time 
steps. In all the cases considered, the proposed methodology produced physically meaningful non- 
negative concentrations. 

4.3. Two-dimensional problem with non-uniform initial condition. This test problem is 
a two-dimensional extension of the problem described earlier in subsection 14.21 The governing 
equations take the following form: 

dc(x,y,t) /<9 2 c(x, y,t) 9 2 c(x,y, t)\ ^ . , n . . _ , . 

+ a 9 = infix := (0,1) X (0,1) X (0,2) (34a) 



dt \ <9x 2 <9y 2 

c(x = 0,y,t) =c(x= l,y,t) =0, c(x, y = 0, t) = c(x, y = 1, t) = (34b) 

J 1 if xe [a, b] x [a, b] 
c(x,y,0) = <^ (34c) 

I otherwise 

Figure [10] gives a pictorial description of the test problem. The analytical solution can be written 
as follows: 

4 1 

c(x, y, t) =— 7T } } (cos(m7ra) — cos(?n7r6))(cos(n7ra) — cos(n7r6)) 

7r z — ' Z — ' mn 

m=l n=l 

sin(n7rx) sin(m7ry) exp [— (m 2 + n 2 )-7r 2 t] (35) 

The numerical convergence of the proposed methodology is shown in Figure [TTJ The performance of 
the proposed methodology with that of the single- field formulation and MATLAB's PDE Toolbox 
is illustrated in Figures [T2] and [T3J. 
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4.4. Two-dimensional problem with anisotropic medium. This problem considers transient 
anisotropic diffusion in a bi-unit square domain. Homogeneous Dirichlet boundary condition is 
applied on the entire boundary. The initial concentration is taken to be zero (i.e., co(x) = 0). The 
diffusivity tensor is taken as follows: 



Four-node quadrilateral element is employed in the numerical simulation. The numerical results 
are generated for two different meshes (XSeed = YSeed = 51 and 101) and two different time steps 
(At = 0.05 s and 0.1 s). Note that XSeed and YSeed denote the number of nodes along x-direction 
and y-direction, respectively. For these cases, the variation of the minimum concentration with time 
under the single- field formulation is shown in Figure [T4"l The proposed methodology produced non- 
negative values for the concentration under all the considered cases, and the minimum concentration 
is zero. In the case of transient isotropic diffusion, the smaller the time step the greater the violation 
of the non-negative constraint [27]- But as evident from Figure [Til this need not be the trend in the 
case of transient anisotropic anisotropic diffusion. Figure [TBI shows the variation of the maximum 
concentration with time under the single-field formulation and the proposed methodology, and there 
is not much difference between the numerical results obtained using the single-field formulation and 
the proposed methodology. Figure [16] compares the contours of the concentration obtained using 
the single-field formulation and the proposed methodology for XSeed = 101 and At = 0.05 s. 
Even for a problem involving anisotropic diffusion, the proposed methodology did not violate the 
non- negative constraint. 



We have presented a novel methodology for transient anisotropic diffusion equations that satisfies 
maximum principles and the non-negative constraint on computational grids with no additional 
restrictions on the time step. The methodology has been developed using the method of horizontal 
lines, and techniques from convex programming. We have shown that the semi-discrete procedure 
based on the standard single-field formulation gives unphysical negative concentrations and violates 
maximum principles. Using several representative numerical examples we have shown that the 
proposed methodology satisfies maximum principles and the non-negative constraint on general 
computational grids with anisotropic and heterogeneous diffusion. The proposed methodology 
performs gives physically meaningful non-negative concentrations even on coarse compuational 




(36) 



with e = 0.05. The volumetric source is taken as follows: 




(37) 



5. CONCLUDING REMARKS 
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grids and for small time steps. We shall conclude the paper by discussing two possible future 
research endeavors in the area of discrete maximum principles. We also briefly outline potential 
challenges one may have to overcome in addressing these research problems. 

(i) A possible future work is to incorporate advection in addition to diffusion, and devise a non- 
negative methodology for both steady-state and transient advection-diffusion equation. How- 
ever, one cannot directly implement the procedure presented in this paper and in references 
[ISJ 03] for advection-diffusion equation, as the advection term makes the spatial differential 
operator non-self-adjoint. 

(ii) Another interesting research problem is to devise a non-negative methodology for both steady 
and transient nonlinear diffusion-type equations. The obvious challenges will be handling 
nonlinearity, and to ensure that the computational cost in obtaining non-negative solutions is 
not prohibitively expensive. 



We now discuss other possible ways of implementing the methods of horizontal and vertical 
lines for transient diffusion- type equations. We will also provide reasons why these approaches 
may not satisfy maximum principles and the non-negative constraint. This discussion will shed 
light on the rationale behind the proposed methodology, and can guide future efforts in developing 
robust solvers for other important parabolic partial differential equations (e.g., transient diffusive- 
reactive systems). All the approaches presented in this appendix employ trapezoidal family of time 
integrators, which can be written as follows: 



where 7 G [0, 1] . (Recall that the parameter 7 used in Section [3] is different from the parameter in 
trapezoidal family of time integrators.) The discussion and conclusions in this appendix will hinge 
on the following result from Matrix Algebra. Given any vector 6^0, the solution of a system of 
linear equations of the form 



will be non- negative (i.e., x y 0) if and only if the matrix A is a monotone. (Recall that ^ denotes 
component- wise inequality.) A matrix is called a monotone if the matrix is invertible and all the 
entries of its inverse are non-negative. For further details on monotone matrices refer to the classic 



6. APPENDIX 




(38) 



Ax = b 



(39) 



texts [22 El EH]. 



is 



6.1. Method of vertical lines at integral time steps. In this paper, this method is referred 
to as the standard single-field formulation. This is the most commonly used method for solving 
transient diffusion equation, and can be found in many introductory texts on finite element methods 
(e.g., |3U| [q"S , 62J). The method is based on standard semi-discrete methodology and Galerkin 
formalism. The corresponding weak form reads: Find c(x, t) G Vt such that we have 

/ w(x ) dc(x,t) dn+ I g^^] . D (x)grad[c(x,t)] dft 
Jn Jn 

= [ w(x) f(x,t) dn+ [ to(x) g p (x, t) dr Vw(x) G Q (40) 

where 

Pi := {c(x, t) G tf 1 (0) | c(x, t) = cp(x, i) on T D } (41) 

and the function space Q is defined previously in equation (|21b|) . After spatial discretization using 
the finite element method, one obtains a system of ordinary differential equations of following form: 

C^L + Kc(t) = f(t) (42) 

The capacity matrix C is symmetric and positive definite, and all the entries of the matrix are non- 
negative. The matrix K is symmetric and positive semi-definite. More importantly, the matrix K 
will not be a monotone if the medium (i.e., the diffusion process) is not isotropic. (If the medium 
is isotropic, it is easy to check that the matrix K is diagonally dominant, and hence it will be a 
monotone matrix.) If a time stepping scheme from the trapezoidal family is employed to solve the 
above ordinary differential equations, one can obtain a system of linear equations of the following 
form: 

{^Kt C + K ) ° {n+1) = f(n+1) + ^Kt C ( cW + ~ 7 ^ (n) ) ( 43 ) 

There are two potential scenarios that can contribute to the violation of the non-negative con- 
straint and maximum principle under the method of vertical lines at integral time steps. Firstly, 
the vector on the right side of equation (|43|) need not be non-negative, as there is no physical con- 
straint requiring that u (n) should be non-negative. Even if the volumetric source is non-negative 
(i.e., /( n+1 ) y 0), c^") y 0, 7 > 0, At > 0, and all the entries of the capacity matrix are non- 
negative; the resulting vector on the right side of the above equation need not be non-negative. One 
possible exception is when 7 = 1 (that is, when the backward Euler is employed). Secondly, the 
matrix on the left side of equation (|43|) may not be a monotone. Even for an isotropic medium, the 
matrix will be monotone only if the time step is greater than a critical time step or by employing 

lumped capacity matrix. Based on the above discussion, the sufficient conditions for the method of 
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vertical lines at integral time levels to satisfy maximum principles and the non-negative constraint 
are as follows: 

• isotropic diffusion, 

• low-order finite elements, 

• backward Euler scheme (i.e., 7 = 1), 

• lumped capacity matrix, 

• select a time step greater than the critical time step, and 

• place constraints on the mesh and element shapes (e.g., well-centered triangular elements, 
rectangular elements with aspect ratio between l/v2 and \/2). 

It is important to note that the above conditions are too restrictive to be able to obtain physically 
meaningful results for practical problems. But this method is commonly employed in many numer- 
ical simulations, and in many commercial finite element packages. Few other remarks about this 
method are in order. 

Remark 6.1. For a discussion on necessary constraints on a finite element mesh to satisfy max- 
imum principles and the non-negative constraint, see references \16\ [14] [29] l45| [43] . However, all 
these constraints are for isotropic diffusion. It is noteworthy that, in the case of anisotropy, a 
computational mesh may not even exist that will ensure the satisfaction of maximum principles and 
the non-negative constraint. 

Remark 6.2. Several studies derived critical time steps with respect to maximum principles. For 
example, see references |60[ [32] . But these derivations for critical time steps are restricted to one- 
dimensional problems, isotropic diffusion, and backward Euler. 

Remark 6.3. It is noteworthy that there is no obvious way of modifying the non-negative formula- 
tions that has been shown recently shown to be successful for steady- state diffusion equations (e.g., 
see references |45} 143] ) to obtain a non-negative formulation for transient diffusion equation under 
the method of vertical lines at integral time steps. This is the reason why this method has not been 
considered as the basis in Section^ 

6.2. Method of horizontal lines at integral time steps. By applying the method of horizontal 

lines at integral time levels and eliminating ?/ n+1 )(x) using the time discretization of trapezoidal 
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family given by equation (|38p . the time discretized equations take the following form: 

^c (n+1) (x) - div[D(x)grad[c (n+1) ]] = / (n+1) (x) + (c (n) (x) + (1 - 7 )Ati/ n )(x)) in U 

(44a) 

c ( n+1 )(x) = 4 n+1 )(x) onr D (44b) 
n(x) • D(x)grad[c (n+1) (x)] = (x) on T N (44c) 

In going from equations (|3ap - (|3d|) to equations (|44ap - (|44cp . the temporal discretization may not 
preserve the non-negative constraint, which should be interpreted in the following sense. One may 
not get a non-negative solution under equations (|44ap - (|44cp even when the solution to the original 
time continuous problem given by equations (|3aj) (|3d|) is non-negative. This is again due to the fact 
that the right side of equation (|44ap can be negative, as there is no physical constraint requiring 
that the rate of concentration v^ n \x) should be non-negative. However, it does not mean that 
the time discrete equation does not satisfy maximum principles and the non-negative equation. 
The above equation is diffusion with decay, and as mentioned earlier, this equation also satisfies 
maximum principles and the non-negative constraint. But, the requirement for the non- negative 
constraint is that / (n+1) (x) + (c( n )(x) + (1 - -f)Atv^ n \x)) > 0. 

6.3. Method of horizontal lines at weighted time levels. We shall perform temporal dis- 
cretization at the weighted time level i n + 7 , which gives rise to the following equations: 

1 - C ( n+7 ) (x) - div[D(x)grad[c( n+7 )]] = f^ n+ ^ (x) + — r (x) in ft (45a) 



jAt jAt 
c< n+ T>(x) = Cp (x,t n+7 ) onr D (45b) 

n(x) • D(x)grad[c( n+7 )] = 4" +7) ( x ) on pN ( 4 5c) 

One can obtain nodal concentrations at weighted time levels (i.e., c( n+7 )) by employing the optimization- 
based solver presented in Section [3J Noting the results presented in Reference [H] on stability issues 
associated with numerical time integration of differential- algebraic equations, the concentration at 
integral time levels is approximated in terms of corresponding quantities at weighted time levels. 
The interpolation scheme is pictorially described in Figure El and can be mathematically written 
as follows: 

c (n+l) = 7C (n+ 7 ) + (! _ 7 ) c («+l+7) ( 46 ) 

The rate of concentration at weighted time levels can be calculated as follows: 

At v ; 
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Figure 2. The left part of the figure shows the usual way of interpolating quanti- 
ties at integral time levels to obtain the corresponding quantities at weighted time 
levels. That is, c^ n+ "^ = (1 — 7)0^"^ + 7c( n+1 ). The right part of the figure shows 
the interpolation of quantities at weighted time levels to obtain the correspond- 
ing quantities at integral time levels, which is adopted in this paper. That is, 
c ( n ) = ■j C ( n - 1 +"/) -f (1 — -y)c( n+1 \ The interpolated quantities are indicated using 
hollow circles. 



The corresponding quantity at integral time levels are calculated as follows: 

v (n+l) = 7U (n+ 7 ) + _ 7 yn+l+ 7 ) ( 48 ) 

The interpolation given by equation (|46p is different from the usual way of interpolating the quan- 
tities at weighted time levels in terms of integral time levels. That is, 

c (n+ 7 ) = ^ _ 7 ) c (n) + 7C (n+l) ( 49 ) 

Figure [2] compares both these interpolation schemes. The only drawback of the method presented 
in this subsection is that it is not self-starting, as we do not have c( n-1+7 ) when n = 1 unless 7 = 1. 
But this drawback can be easily overcome by employing the backward Euler scheme (i.e., 7 = 1) 
for the first time level, and then employ the method for subsequent time levels. Therefore, the 
method presented in this subsection can be considered as an alternate to the method presented in 
Section [3] to satisfy maximum principles and the non- negative constraint for transient diffusion- type 
equations. 
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FIGURE 3. One dimensional problem with uniform initial condition: This figure 
shows the concentration at x = 0.6 as a function of time for 7 = 1 (top figure) and 
7 = 0.1 (bottom figure). The time step is taken as At = 0.001 s, and five equally 
spaced linear finite elements are employed. The numerical solutions obtained from 
the single-field formulation and the proposed methodology are compared with the 
analytical solution. From the maximum principle, it is known that the analytical 
solution is bounded above by unity. The numerical solution from the single-field 
formulation exceeds unity while the proposed methodology satisfies the maximum 
principle. 
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Figure 4. One-dimensional problem with non- uniform initial condition: This fig- 
ure compares the concentration obtained using the proposed methodology with the 
analytical solution at various instants of time. For this test problem, the solution 
should be between zero and unity. The time step used in the numerical simulation is 
At = 10~ 4 s. As one can see from the figure, the proposed methodology performed 
well, and it did not violate the maximum principle and the non- negative constraint. 
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versity, College Station, Texas 77843. 
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-log(h) -log(h) 
(a) single-field formulation (b) proposed methodology 

Figure 5. One-dimensional problem with non-uniform initial condition: This figure 
compares the numerical convergence of the single-field formulation and the proposed 
methodology with simultaneous spatial and temporal refinements such that At oc 
(Ax) 2 . In this numerical simulation, we have taken 7=1, and X = 0.5 s. The 
coarsest mesh has 100 nodes, and the corresponding time step used for this mesh is 
At = 0.05 s. The terminal rates of convergence in Z/2-norm and i7 1 -seminorm are 
also shown in the figure. 
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FIGURE 6. One-dimensional problem with non-uniform initial condition: This figure 
shows the variation of the minimum concentration with time for different meshes 
under the single-field formulation. The time step is takes to be At = 0.001 s. The 
mesh is discretized using equally spaced nodes. The single-field formulation produces 
unphysical negative values for the concentration. However, for this one-dimensional 
problem, the extent of the violation of the non-negative constraint decreases with 
the mesh refinement. (It should be noted that the violation may not decrease with 
mesh refinement if the diffusion is anisotropic 
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Figure 7. One-dimensional problem with non-uniform initial condition: This figure 
shows the variation of the maximum concentration with time for different meshes 
under the single-field formulation. The time step is takes to be At = 0.001 s. The 
mesh is discretized using equally spaced nodes. The single-field formulation violates 
the maximum principle. However, for this one-dimensional problem, the extent of 
the violation of the maximum principle decreases with the mesh refinement. (It 
should be noted that the violation may not decrease with mesh refinement if the 
diffusion is anisotropic 
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Figure 8. One-dimensional problem with non- uniform initial condition: This fig- 
ures illustrates that, for a given mesh, the extent of the violation of the maximum 
principle and the non-negative constraint increases with a decrease in the size of 
the time step under the single-field formulation. The top figures are obtained us- 



ing At = 10 4 s and the bottom figures are obtained using At 



10' 



s. 



The 



left figures are obtained using a computational mesh of 25 equally spaced nodes, 
and the right figures are obtained using a computational mesh of 42 equally spaced 
nodes. The time of interest is taken as I = 10 s. For the mesh with 25 nodes, 
the percentage of nodes violated the maximum principle is 8% for both At = 10~ 4 s 
and At = 10 -7 s; and the percentage of nodes violated the non-negative constraint 
is 32% for At = 10" 4 s and 40% for At = 10" 7 s. For the mesh with 42 nodes, 



the percentage of nodes violated the maximum principle is 0% for At 



10' 



and 9.52% for At = 10 s; and the percentage of nodes violated the non-negative 
constraint is 0% for At = 10~ 4 s and 38.1% for At = 10 -7 s. 
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Figure 9. One-dimensional problem with non-uniform initial condition: This figure 
illustrates the effect of mesh refinement and small time steps on the performance of 
the proposed methodology. The top figures are obtained using At = 10~ 4 s and the 
bottom figures are obtained using At = 10~ 7 s. The left figures are obtained using 
a computational mesh of 25 equally spaced nodes, and the right figures are obtained 
using a computational mesh of 42 equally spaced nodes. The time of interest is 
taken as X = 10~ 4 s. The proposed methodology performed well even on coarse 
meshes, and for small time steps. 



(0,1) 



c(x,y = 1,0) = 



o 

II 



(0,0) 





___(&,&) 


(a, a) 


(6, a) 



(1,1) 



o 



;i,o) 



c(x,y = 0,0) = 

Figure 10. Two-dimensional problem with non- uniform initial condition: A picto- 
rial description of the problem described in subsection !4.31 The shaded region has an 
initial concentration of c(x, y, t = 0) = 1, and the remaining part of the domain has 
an initial condition of c(x, y,t = 0) = 0. Homogeneous Dirichlet boundary condition 
is prescribed on the entire boundary. 
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Figure 1 1 . Two-dimensional problem with non-uniform initial condition: This fig- 
ure illustrates the numerical convergence of the proposed methodology. We have 
taken 7 = 1, and the length of the time interval is taken as I = 0.3 s. A hierar- 
chy of meshes are employed in the numerical study. The initial mesh has 31 nodes 
along each direction, and the initial time step is taken as At = 0.01 s. The mesh 
and the time step are simultaneously refined as At oc (Ax) 2 . The terminal rates of 
convergence in L2-norm and i^ 1 -seminorm for the proposed methodology are 1.50 
and 1.49, respectively. 
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(c) Violation of the maximum principle (d) Proposed methodology 

Figure 12. Two-dimensional problem with non- uniform initial condition: This fig- 
ure compares the concentrations obtained from the single-field formulation and the 
proposed methodology with the analytical solution at time level t = At = 10 -4 s. 
Subfigure (b) shows that the single-field formulation violates the non-negative con- 
straint, as 36% of nodes have negative concentrations. The obtained minimum 
concentration is —0.01221. Subfigure (c) shows that the single-field formulation vi- 
olates the maximum principle, as 1% of nodes having concentrations greater than 
unity. The obtained maximum concentration is 1.02039. Subfigure (d) shows that 
the concentration obtained from the proposed methodology satisfies the maximum 
principle, and the non- negative constrain^ In subfigure (b), the regions with nega- 
tive concentrations are indicated in white color. In subfigure (c), the regions with 
concentrations greater than unity are indicated in white color. 
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Figure 13. Two-dimensional problem with non-uniform initial condition: This fig- 
ure compares the numerical solutions from MATLAB's PDE Toolbox and the pro- 
posed methodology at time level t = At = 10 -4 s. Subfigure (a) shows the compu- 
tational mesh used in the numerical simulation. Subfigure (b) shows that numerical 
solution from the MATLAB's PDE Toolbox violates the non- negative constraint, as 
40% of the nodes have negative concentrations. The regions with negative concentra- 
tions are indicated in white color. The obtained minimum concentration is —0.0339. 
Subfigure (c) shows that the numerical solution from MATLAB's PDE Toolbox vi- 
olates the maximum principle, as 1.2% of nodes have concentrations greater than 
unity. The regions with concentrations greater than unity are indicated in white 
color. The obtained maximum concentration is 1.0397. Subfigure (d) shows that 
the proposed methodology satisfies the maximum principle and the non-negative 
constraint on the computational mesh generated by MATLAB. 
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Figure 14. Two-dimensional problem with anisotropic medium: This figures shows 
the variation of the minimum concentration under the single-field formulation. The 
results are shown for two different meshes (XSeed = YSeed = 51 and 101), and 
for two time steps (At = 0.05 s and 0.1 s). Note that XSeed and YSeed denote 
the number of nodes along x-direction and y-direction, respectively. The single-field 
formulation produced negative concentrations for both the meshes and for both the 
time steps. The proposed methodology produced non-negative solutions under all 
the cases considered. 
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FIGURE 15. Two-dimensional problem with anisotropic medium: This figures shows 
the variation of the maximum concentration under the single-field formulation (top 
figure) and the proposed methodology (bottom figure). The results are shown for 
two different meshes (XSeed = YSeed = 51 and 101), and for two time steps (At = 
0.05 s and 0.1 s). Note that XSeed and YSeed denote the number of nodes along 
x-direction and y-direction, respectively. As evident from the figure, the single- 
field formulation and the proposed methodology produced similar results for the 
maximum concentration with respect to time. 
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FIGURE 16. Two-dimensional problem with anisotropic medium: This figure shows 
the contours of the concentration under the single-field formulation (left) and the 
proposed methodology (right) at time = 1 s (top) and time = 2 s (bottom). The 
time step is taken as At = 0.05 s, and XSeed = YSeed = 101. The number of nodes 
along x-direction and y-direction are, respectively, denoted by XSeed and YSeed. 
The regions that have violated the non-negative constraint are indicated in white 
color. 
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